EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH 



PH-EP/2005-057 
25th November 2005 



Measurement of the Strong Coupling as from 
Four-Jet Observables in e + e _ Annihilation 



The OPAL Collaboration 



Abstract 

Data from e + e~ annihilation into hadrons at centre-of-mass energies between 
91 GeV and 209 GeV collected with the OPAL detector at LEP, are used to study 
the four-jet rate as a function of the Durham algorithm resolution parameter y cu t- 
The four-jet rate is compared to next-to- leading order calculations that include the 
resummation of large logarithms. The strong coupling measured from the four-jet 
rate is 

a s (M z o) = 

0.1182 ± 0.0003(stat.) ± 0.0015(exp.) ± 0.0011(had.) ± 0.0012(scale) ± 0.0013(mass), 

in agreement with the world average. Next-to-leading order fits to the D-parameter 
and thrust minor event-shape observables are also performed for the first time. We 
find consistent results, but with significantly larger theoretical uncertainties. 
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1 Introduction 

The annihilation of electrons and positrons into hadrons allows a precise test of Quantum 
Chromodynamics (QCD). Many observables have been devised which provide a convenient 
way of characterizing the main features of such events. Multi-jet rates are predicted in 
perturbation theory as functions of the jet-resolution parameter, with one free parameter, 
the strong coupling as- Events with four quarks in the final state, qqqq, or two quarks and 
two gluons, qqgg, may lead to events with four-jet structure. In leading order perturbation 
theory, the rate of four-jet events in e + e~ annihilation is predicted to be proportional to 
ctg. Thus, the strong coupling is measured by determining the four-jet production rate in 
hadronic events and fitting the theoretical prediction to the data. 

Calculations beyond leading order are made possible by theoretical progress achieved 
during the last few years. For multi-jet rates as well as numerous event-shape distributions 
with perturbative expansions starting at 0(as), matched next-to-leading order (NLO) and 
next-to-leading logarithmic approximations (NLLA) satisfactorily describe data over large 
kinematically allowed regions at many centre-of-mass energies [1-4]. 

Measuring the four-jet rate is not the only way to determine the strong coupling; 
various event-shape observables sensitive to the presence of four partons in the final state 
may also be studied. In addition to measuring the four-jet rates, we also studied two 
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event-shape observables which are sensitive to the four jet structure of the final state, 
the D-parameter [5] and the thrust minor (T min ) [6]. For these observables the NLLA 
predictions are only available in the nearly planar limit. Furthermore, the NLO corrections 
are rather large, which indicates that the neglected higher order terms in the perturbative 
prediction are important [7, 8] and lead to a large uncertainty in the determination of the 
strong coupling. 

In this analysis we use data collected by OPAL at LEP at 13 centre-of-mass energies 
covering the full LEP energy range of 91—209 GeV. A previous publication by the OPAL 
Collaboration performed a simultaneous measurement of «s and the QCD colour factors 
with data taken at 91 GeV [9]. The same theoretical predictions for the four-jet rate are 
used for this analysis, with the colour factors Ca and C-p set to the values expected from 
the QCD SU(3) symmetry group, Ca = 3 and CV = 4/3. More importantly, our analysis 
extends the energy range up to 209 GeV and we also investigate the D-parameter and 
Tmin event-shape observables. Similar results from other LEP collaborations can be found 
in [10,11]. 

The OPAL Collaboration has also published a comprehensive measurement of jet-rates 
in the energy range between 91 and 209 GeV [12] and the distributions of the D-parameter 
and T min amongst many other event shape observables [13]. The data used in the present 
paper are identical to those presented in [12, 13]. However, the determinations of as 
in [12,13] were based on the two-jet rate and event shapes sensitive to three-jet processes; 
in contrast the present paper focuses on different calculations for observables sensitive to 
four-jet production. 



2 Observables 

2.1 The Four- Jet Rate 

Jet algorithms are applied to cluster the large number of particles of a hadronic event into 
a small number of jets, reflecting the parton structure of the event. For this analysis we 
use the Durham scheme [2]. Defining each particle initially to be a proto-jet, a resolution 
variable is calculated for each pair of proto-jets i and j: 

Vv = W2 t 1 ~ cos %)> (!) 

where Ei and Ej are the energies of jets i and j, cos#y is the angle between them and 
E vis is the sum of the energies of the visible particles in the event (or the partons in a 
theoretical calculation). If the smallest value of yij is less than a predefined value y cut , the 
pair is replaced by a new proto-jet with four-momentum p% = pf + pf, and the clustering 
starts again. Clustering ends when the smallest value of y^ is larger than y cut , and the 
remaining proto-jets are accepted and counted as final selected jets. 

The NLO calculation predicts the four-jet rate _R 4 , which is the fraction of four-jet 
events, as a function of as and y cut . The NLO prediction can be written in the following 
way [14]: 

i? 4 (y C ut) = a4 " Jct(ycUt) = V 2 B,(y C ut) + V 3 [C,(y C ut) + 3/2(/3 \ogx, - 1) B 4 (y cut )\ (2) 

Ctot 
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where o" 4 _j et is the cross-section for the production of a hadronic event with four jets at 
fixed Ucut, Ctot the total hadronic cross-section, 77 = as Cf/2tt, x m = /// ' y/s with fi being 
the renormalization scale, \fs the centre-of-mass energy, and /3rj — (11 — 2nf/3) with nf 
the number of active flavours 1 . The coefficients B4 and C4 are obtained by integrating 
the massless matrix elements for e + e~ annihilation into four- and five-parton final states, 
calculated by the program Debrecen 2.0 [14] 2 . Eq. (2) is used to predict the four-jet rate as 
a function of y cnt . The fixed-order perturbative prediction is not reliable for small values 
of ?/ cu t due to terms of the form a§ ln m (y cnt ) which enhance the higher order corrections. 
An all-order resummation of such terms (NLLA) is possible for the Durham clustering 
algorithm [2]. The NLLA calculation is combined with the NLO-prediction using the 
so-called modified R-matching scheme [9]. In this scheme the terms proportional to rf 
and rf are removed from the NLLA prediction _R NLLA and the difference is then added to 
the NLO calculation i? 4 from Eq. 2 with the result: 

^R-match = ^NLLA + _ ^NLLA) + _ ^NLLA _ 3 / 2 (B A _ ^NLLA)^ (3) 

where _g NLLA and (7 NLLA a re the coefficients of the expansion of i? NLLA , taking the same 
form as Eq. (2), and the term and y cut dependence have been suppressed for clarity. 



2.2 Thrust minor (T m [ n ) and D-parameter 

The properties of hadronic events may also be described by event-shape observables. 
These are used to characterize the distribution of particles in an event as "pencil-like", 
planar, spherical, etc. The variable T min and the D-parameter are sensitive to hadronic 
configurations with more than three partons in the final state. These observables are 
defined as follows: 

Thrust minor T min : The thrust is defined by the expression [6] 

T = max — , (4) 



n 



Ei \p% 



where pi is the three-momentum of particle % and the summation runs over all 
particles. The thrust axis Ht is the direction n which maximizes the expression 
in parentheses. In a further step the maximization in Eq. (4) is performed with 
all momenta projected onto the plane perpendicular to ny. The resulting vector is 
called ^T ma or - At the last step the expression in parentheses in Eq. (4) is evaluated 
again for the vector nT min which is perpendicular to both Ht and riT ma . or , this is T n 



miri' 



£)-parameter: This observable is based on the linearized momentum tensor 

n»/J ^2i(P?Pi)/\Pi\ y~. Q 1 O Q (K\ 

o = ^ | | , a, /J =1,2,3, (5) 

Ei\Pi\ 

1 In this analysis the number of active flavours is set to five. 

2 The Durham y cut values are chosen to vary in the range between 0.00001 and 0.3162, similar to the 
study in Ref. [14]. 



6 



where the sum runs over particles, i, and a and (3 denote the Cartesian coordinates 
of the momentum vectors. The three eigenvalues Xj of this tensor define D through 



D = 27AiA 2 A 3 . 



(6) 



The NLO calculation predicts the event-shape observable distributions 



1 da 



V 2 Bo 4 (Oa) + V 3 [Co 4 (Oa) + 3/2(/3 log^ - 1) B 0a {O a )] 



(7) 



(T t ot dOi 



with 4 being the D-parameter or T min . The coefficients B 0a and Co A are obtained 
by integrating the massless matrix elements for e + e~ annihilation into four-parton final 
states, calculated by the program Debrecen 2.0. 



Although the data used in the present paper are identical to those presented in [12,13], for 
the sake of completeness, and in order to present a coherent discussion of the systematic 
uncertainties, we give a brief account of the analysis procedure below. 



A detailed description of the OPAL detector can be found in Ref. [15]. This analysis 
relies on the reconstruction of charged particle trajectories and on the measurement of 
energy depositions in the electromagnetic calorimeters. Tracking of charged particles was 
performed with the central detector, which was located in a solenoidal magnet providing 
an axial magnetic field of 0.435 T. The central detector contained a silicon microvertex 
detector, and three drift chamber devices: an inner vertex chamber, a large volume jet 
chamber and surrounding z-chambers to measure the z-coordinate 3 . Most of the tracking 
information was obtained from the jet chamber, which provided up to 159 measured 
space points per track, and almost 100% tracking efficiency in the region | cos#| < 0.98. 
Electromagnetic energy was measured by lead glass calorimeters surrounding the magnet 
coil, separated into a barrel (|cos6>| < 0.82) and two end-cap (0.81 < |cos#| < 0.98) 
sections. The electromagnetic calorimeter consisted of 11704 lead glass blocks with a 
depth of 24.6 radiation lengths in the barrel and more than 22 radiation lengths in the 
end-caps. 

3.2 Data Samples 

The data used in this analysis were collected by the OPAL detector between 1995 and 2000 
and correspond to 14.5 pb~ of 91 GeV data, 10.4 pb~ of LEP1.5 data with centre-of- 
mass energies of 130 GeV and 136 GeV and 706.4 pb" 1 of LEP2 data with centre-of-mass 

3 In the OPAL right-handed coordinate system the x-axis pointed towards the centre of the LEP ring, 
the y-axis pointed approximately upwards and the z-axis pointed in the direction of the electron beam. 
The polar angle and the azimuthal angle 4> are defined with respect to z and x, respectively, while r is 
the distance from the z-axis. 



3 Analysis Procedure 



3.1 The OPAL Detector 
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energies ranging from 161 to 209 GeV. The highest energy runs at average centre-of- 
mass energies between 192 and 207 GeV have an energy spread of about 1-2 GeV and 
we grouped these into six main energy points. The data at 91 GeV were taken during 
calibration runs at the Z° peak at the beginning of each year and during data-taking in 
1996-2000. The breakdown of the data samples, the mean centre-of-mass energies, the 
energy range, the data-taking years and the collected luminosities are given in Table 1. For 



average 


energy 


year 


luminosity 


events 


events 


energy in GeV 


range in GeV 




(Pb" 1 ) 


selected 


predicted 


91.3 


91.0-91.5 


1996-2000 


14.7 


395695 


— 


130.1 


129.9-130.2 


1995, 1997 


5.31 


318 


368.4 


136.1 


136.0-136.3 


1995, 1997 


5.95 


312 


329.7 


161.3 


161.2-161.6 


1996 


10.1 


281 


275.3 


172.1 


170.2-172.5 


1996 


10.4 


218 


232.2 


182.7 


180.8-184.2 


1997 


57.7 


1077 


1083.5 


188.6 


188.3-189.1 


1998 


185.2 


3086 


3130.1 


191.6 


191.4-192.1 


1999 


29.5 


514 


473.0 


195.5 


195.4-196.1 


1999 


76.7 


1137 


1161.3 


199.5 


199.1-200.2 


1999, 2000 


79.3 


1090 


1130.8 


201.6 


201.3-202.1 


1999, 2000 


37.8 


519 


526.5 


204.9 


202.5-205.5 


2000 


82.0 


1130 


1089.6 


206.6 


205.5-208.9 


2000 


138.8 


1717 


1804.1 



Table 1: The average center-of-mass energy, the energy range, the year of data taking 
and the integrated luminosity for each data sample, together with the numbers of selected 
data and number of events expected from Monte Carlo simulation. The horizontal lines 
separate the four energy intervals. 

presentation purposes we combine the data above 91 GeV into three energy points [13]. In 
the combination the individual data sets are weighted by the product of their luminosity 
and signal cross-section. The LEP1.5 data samples provide an energy point at 133.3 GeV, 
while the LEP2 samples give energy points at 177.4 GeV and 197.0 GeV corresponding 
to the energy ranges from 161 to 184 GeV and from 188 to 209 GeV. 

3.3 Monte Carlo Samples 

The same samples of Monte Carlo simulated events as in the analysis of [13] are used to 
correct the data for experimental effects such as acceptance, resolution and backgrounds. 
The process e + e _ — > qq is simulated using JETSET 7.4 [16] at v /i=91.2 GeV, and us- 
ing KK2f [17] with fragmentation performed by PYTHIA 6.125 [18] at higher energies. 
Corresponding samples using HERWIG 6.2 [19] are used for systematic checks. Elec- 
troweak four-fermion background processes are simulated using KORALW 1.42 [20] with 
GRC4f [21] matrix elements and fragmentation performed by PYTHIA. The Monte Carlo 
samples generated at each energy point are processed through a full simulation of the 
OPAL detector [22] and reconstructed in the same way as the data. In addition, for 
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comparisons with the corrected data and to estimate fragmentation effects, large sam- 
ples of generator-level Monte Carlo events are employed, using the parton shower models 
PYTHIA 6.158, HERWIG 6.2 and ARIADNE 4.11 [23]. Each of these fragmentation mod- 
els contains a number of tunable parameters; these were adjusted by tuning to previously 
published OPAL data at ^ ~ 91 GeV as described in Ref. [24] for PYTHIA/ JETSET 
and in Ref. [25] for HERWIG and ARIADNE. 



3.4 Selection of Events 

The selection of events for this analysis follows exactly the study described in [13]. It 
consists of three main stages: the identification of hadronic event candidates, the removal 
of events with a large amount of initial state radiation (ISR), and the removal of four- 
fermion background events. 

The selection of hadronic events is based on cuts on particle multiplicity mainly to 
remove leptonic final states and visible energy and longitudinal momentum balance mainly 
to remove two-photon events, e + e~ — > e + e~7*7*, with hadronic final states. 

Standard criteria [13] are used to select well measured tracks and clusters of energy 
deposits in the calorimeter. The number of good charged particle tracks is required to be 
greater than six. The polar angle of the thrust axis is required to satisfy | cos(#t)| < 0.9 
in order that the events be well contained within the detector acceptance. 

The effective e + e~ centre-of-mass energy after excluding all ISR photons, y/~p, is esti- 
mated for each selected event using the algorithm described in Ref. [26]. At centre-of-mass 
energies of 130 GeV and above, we require that y/s — yfF < 10 GeV in order to select 
non-radiative events. 

At energies above the W + W~ production threshold, electroweak four-fermion events, 
especially those involving qqqq final states, become a substantial background. These are 
reduced by using the standard OPAL W + W~ likelihood-based selections [27]. At centre- 
of-mass energies of 161 GeV and above, the W + W~ — > qqqq likelihood is required to 
satisfy £qq qq - < 0.25 and the W + W~ — > qqiu likelihood is required to satisfy £ qq ^ < 0.5. 

The numbers of selected candidate events obtained after applying the selection cuts 
are given in Table 1. The numbers are consistent with expectations based on Monte Carlo 
simulation 4 . After all cuts, the acceptance for non-radiative signal events 5 ranges from 
88.5% at 91 GeV, where the loss in acceptance is mostly due to the cut on the polar 
angle of the thrust axis, to 76.5% at 207 GeV, where additional cuts on y/F and four- 
fermion background rejection reduce the efficiency. The residual four-fermion background 
is negligible below 161 GeV, and increases from 2.1% at 161 GeV to 6.2% at 207 GeV. 



4 Thc numbers of events expected on the basis of Monte Carlo simulations are given in all cases except 
for 91 GeV; to perform an accurate prediction close to the Z° peak would require a much more careful 
investigation of the beam energy and luminosity than is required for the present analysis. 

5 Defined for this purpose as those fulfilling yfs — y/F < 1 GeV, with y/F the effective e + e~ centre-of- 
mass energy after excluding all ISR photons. 
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3.5 Corrections to the data 

In order to mitigate the effects of double counting of energy in tracks and calorimetry, a 
standard algorithm is adopted [13] which associates charged particles with electromagnetic 
calorimeter clusters, and subtracts the estimated contribution of the charged particles 
from the calorimetry energy. All selected tracks, and the electromagnetic calorimeter 
clusters remaining after this procedure, are used in the evaluation of the distributions. The 
resulting distributions with all selection cuts applied are called detector-level distributions. 

The expected number of remaining four-fermion background events, Q, is subtracted 
from the number of data events, iVj, for each data point of each distribution 6 . The effects 
of detector acceptance and resolution and of residual ISR are then accounted for at each 
data point. Since in the present analysis the Monte Carlo model gives a good description 
of the data, and migration between data points is small, this procedure is justified. 

Two distributions are formed from Monte Carlo simulated signal events for each 
observable; the first, at the detector level, treats the Monte Carlo events identically 
to the data, while the second, at the hadron level, is computed using the true four- 
momenta of the stable particles 7 in the event, and is restricted to events whose s' satisfied 
y/s — y/s' < 1 GeV. The Monte Carlo ratio of the hadron level to the detector level for 
each data point, ccj, is used as a correction factor for the data, yielding the corrected data 
point Ni = ai(Ni - Q). 

The hadron-level distribution is then normalized to unity: Pi = Ni/N, where the sum 
N = J2k Nk includes the underflow and overflow data points where appropriate. 

3.6 Fit Procedure 

Our measurement of the strong coupling as is based on fits of QCD predictions to the 
corrected distributions. The theoretical predictions as described in Section 2 provide 
distributions at the parton level. The parton-level distributions are obtained from the 
partons after the parton shower, just before the hadronization. In order to confront the 
theory with the hadron-level data, it is necessary to correct for hadronization effects. 
This is done by calculating the distributions at both the hadron and the parton level 
using PYTHIA and, as a cross-check, with the HERWIG and ARIADNE models. The 
theoretical prediction is then multiplied by the ratio of the hadron- and parton-level 
distributions. 

A x 2 -value for each energy point is calculated using the following formula 

n 

X 2 = E(04,i - 4 (a s )f lc °)(^(0 4 ))- 1 (0 4j - 0« heo ), (8) 

with i,j running over all data points in the fit range and 4 corresponding to either the 
four-jet rate R4, the differential distribution T m i n or the D-parameter, while VijiO^) is the 
covariance matrix. The \ 2 value is minimized with respect to as f° r each energy point 
separately. 

6 A data point is the center of a bin of a distribution for the D-parameter and T m i n , while for the 
four-jet rate we have points at discrete values of y cu t where the analysis is performed. 

7 All charged and neutral particles with a lifetime longer than 3 x 10~ 10 s are treated as stable. 
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The four-jet rate is an integrated distribution, while the D-parameter and T min dis- 
tributions are differential ones. Therefore in the four-jet rate distribution a single event 
can contribute to several y cut -data point and for this reason the data points are corre- 
lated. The complete covariance matrix Wij is determined from four-jet rate distributions 
calculated at the hadron level. Subsamples are built by choosing 1000 events randomly 
out of the set of all generated Monte Carlo events. A single event may be included in 
several subsamples, but the impact on the final covariance matrix is expected to be small 
and, therefore, is neglected [28]. For every energy point 1000 subsamples are built. The 
covariance matrix is then used to determine the correlation matrix, = Wij/di<7j, with 
&i — y/Wu. The covariance matrix V^R^) used in the x 2 fit is then determined using 
the statistical error <7j of the data sample at data point i and the correlation matrix 

Pij '■ = Pij^iVj- 

For the fit to the event-shape observables the covariance matrix Vij{Oi) for Pj is 
computed by transforming the diagonal Poisson covariance matrix for the uncorrected 
data points JVj [13]: 

v* = £ ^§r Nk = w^ alNk ( mk ~ % ) ( N5jk ~ ^) • (9) 

k k ^ 



3.7 Combination of Energy Points 

For presentation purposes the data of several energy points are combined. The four-jet 
rate and event shape distributions measured at the different energy points are averaged 
using the products of luminosity and cross-section. The fit results for a s , however, are 
combined using the procedure of Ref. [13]. 

In brief the method is as follows. The «s measurements to be combined are first 
evolved to the common scale of the combination, Q = y/so, assuming the validity of 
QCD. The measurements are then combined using a weighted mean method, based on 
minimizing the x 2 between the combined values and the measurements. If the measured 
values evolved to a common scale Qo = ^/so are denoted with covariance matrix V, 
the combined values, «s(v^o); are gi ven by 

«s(v / so) = H w ^ where w i = JHV// -i\ > ( 10 ) 

i Lj,kK v )jk 

and i,j enumerate the individual measurements. Only experimental systematic errors 
(assumed to be fully correlated, V-j = o^a'^ between measurements) are taken to con- 
tribute to the off-diagonal elements in V. All error contributions (statistical, experimen- 
tal, hadronization, scale uncertainty and the uncertainty from massless calculations) are 
taken to contribute to the diagonal elements. The hadronization and scale uncertain- 
ties are computed from the as values obtained with the alternative hadronization models 
and with the upper and lower values of the renormalization scale, respectively, and com- 
bined according to Eq. 10. The uncertainty for the massless calculations is also combined 
according to Eq. 10. 
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4 Systematic Uncertainties 



Several sources of possible systematic uncertainty are studied as in [13]. All systematic 
uncertainties are taken as symmetric. 

Experimental uncertainties: Contributions to the experimental uncertainties are 
estimated by repeating the analysis with varied cuts or procedures. For each systematic 
variation the value of as is determined and then compared to the result of the standard 
analysis (default value). In each case, the difference with respect to the default value is 
taken as a systematic uncertainty. 

(1) The algorithm to avoid the double counting of energy in tracks and calorimetry 
is not applied. Instead all observables are computed using all charged tracks and 
electromagnetic calorimeter clusters. The effects of double counting are then taken 
into account exclusively through the detector correction. 

(2) The containment cut is tightened to be well within the barrel, | cos(#t)| < 0.7. 

(3) Instead of using PYTHIA for the correction of detector effects as described in Sec- 
tion 3.5, events generated with HERWIG are used. 

(4) The fit range 8 is changed. Two different cases are considered. First the fit range 
is decreased by one data point at each edge of the fit range. Second the fit range 
is extended by one data point at each edge of the fit range. The larger deviation 
from the default fit is taken as a systematic uncertainty. In order to take statistical 
fluctuations into account, the deviation is calculated using the average deviation of 
a fit applied to 50 Monte Carlo samples. 

(5) The algorithm to compute s' is replaced by a simpler version accounting for at most 
one initial-state photon [26]. 

(6) The cut on the likelihood variable £ q qqq to reject background from W + W~ — > qqqq 
events is changed from 0.25 to 0.1 and 0.4. The larger deviation from the default 
value is taken as the systematic uncertainty. 

(7) The cut on the likelihood variable C q q£ u to reject background from W + W~ — > qqiu 
events is changed from 0.5 to 0.25 and 0.75. The larger deviation from the default 
value is taken as the systematic uncertainty. 

(8) The amount of subtracted four-fermion background is varied by ±5%. The larger 
deviation from the default value is taken as the systematic uncertainty. 

Variation (5) applies only to the data taken above the Z° resonance and variations (6)- (8) 
only to data taken at and above the W-pair threshold at 161 GeV. All experimental uncer- 
tainties are added in quadrature and the result is quoted as the experimental systematic 
uncertainty. None of the variations contribute dominantly to the overall experimental 
systematic uncertainty. 

8 The determination of the fit range is explained in section 5.1.1 
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Hadronization: The uncertainties associated with the hadronization correction are 
evaluated by using HERWIG and ARIADNE instead of PYTHIA. The larger change in 
as resulting from these alternatives is taken as the error. 

Scale uncertainties: The uncertainty, associated with missing higher order terms in 
the theoretical prediction, is assessed by varying the renormalization scale factor x M . The 
predictions of a complete QCD calculation would be independent of x^, but a finite-order 
calculation such as that used here retains some dependence on x^. The renormalization 
scale x^ is set to half and twice the value of x^ from the default fit. The larger deviation 
from the default value of as is taken as the systematic uncertainty. The uncertainty band 
method as discussed in [29] cannot be applied to this analysis, because the rescaling of 
the NLLA terms and the use of different matching schemes are currently not possible. 
However, the renormalization scale variation in the range 0.5 < x^ < 2.0 gives scale uncer- 
tainties for most observables consistent with uncertainties obtained with the uncertainty 
band method [29]. 

Uncertainties from massless QCD calculations: In the case of the four-jet rate an 
additional uncertainty from massless QCD calculations is considered. The QCD prediction 
used is only available for the massless case. Effects from massive b quarks are expected and 
the corresponding uncertainties on the strong coupling are evaluated for the four-jet rate. 
The overall four-jet rate can be written as the sum of the the four-jet rate originating from 
b quark and from light quark events, i? 4 = fbR\ + (1 — /b)-R 4 ght = (/b-B 4 + (1 — /b))-R 4 ght , 
with /b being the fraction of b quark events, R\ the four-jet rate originating from b quark 
events, i? 4 ght the four-jet rate originating from light quark events and 5 4 = R\/R l l &ht . 
The Standard Model expectation for the fraction /b computed in [30] varies between 
0.216 for data taken at LEP1 and 0.169 at 207 GeV. Following the studies of QCD heavy 
mass effects [31] B 4 ps 0.9 at y/s = M z o, leading to a change in the overall four-jet rate 
i? 4 of « 2%. Even though the effect is expected to decrease with increasing y/s, a value 
of B4 = 0.9 is used as a conservative estimate for the whole energy range. Since i? 4 ~ 
in LO QCD we have Aas/^s — Ai? 4 /2i? 4 and thus we expect a ~ 1% relative change in 
as- This is taken to be our massless QCD uncertainty. 

5 Results 

5.1 Four- Jet Rate Distributions 

The four-jet rates for the four energy intervals, after subtraction of background and cor- 
rection for detector effects, are shown in Figure 1. Superimposed we show the distribution 
predicted by the PYTHIA, HERWIG and ARIADNE Monte Carlo models. In order to 
allow for a more quantitative comparison between data and models, the inserts in the 
upper right corner show the differences between data and each model, divided by the 
combined statistical and experimental error in each data point The sum of squares of 
these differences would, in the absence of correlations, represent a x 2 between data and 
the model. However, since correlations are present, such x 2 values should be regarded 
only as a rough indication of the agreement between data and the models. The three 
models are seen to describe the high energy data well. Some discrepancies are, however, 
seen in the 91 GeV data which have a much better precision. The different Monte Carlo 
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models do not agree with each other; ARIADNE tends to give the best description of the 
data. 

5.1.1 Determination of ct s 

For the QCD prediction of Ri{y C ut) the combined C(o;g)+NLLA calculation (Eq. 3) is 
used. The scale parameter is set to one. The fit range is 0.0028 < y cut < 0.0371 for 
data taken at 91 GeV, 0.0018 < y cut < 0.0178 for data taken at LEP1.5 and 0.0018 < 
Vent < 0.0056 for data taken at LEP2. The fit range is determined by requiring that the 
hadronization corrections be less than 20% and the detector corrections be less than 50%. 
As shown in [14] the theoretical uncertainty due to higher order missing terms increases 
significantly for y cut values below 0.005. Therefore we impose a further constraint on 
the y cut region - that the theoretical prediction should not vary by more than 5% as the 
renormalization scale factor x M is varied from 0.5 to 2 9 . The fit range for the LEP2 
data contains only three data points; an enlarged fit range would lead to only a minor 
gain in statistical precision. We require the same fit range for all energy points above 
136 GeV. In Figure 2 the hadron-level four-jet rate distributions for the four energy 
intervals are shown together with the fit result. The fit ranges cover the falling slope 
of the distributions at large y cut , where the perturbative QCD predictions adequately 
describe the data. The increasing slope of the distributions at low y cnt values is less 
well described by the perturbative QCD description. In this region the hadronization 
corrections become large (more than 100%) and the various hadronization models start 
to differ significantly. The fits performed at individual energy points are combined in 
the four different energy bins: 91.2 GeV; 130 and 136 GeV; 161, 172 and 183 GeV; and 
189-209 GeV using the method described in Section 3.7. The results are summarized in 
Table 2 and shown in Figure 3. At LEP1 the statistical uncertainty is much smaller than 
the theoretical one. The hadronization scales as an inverse power of y/s and therefore the 
hadronization uncertainty decreases with increasing energy. The fit range for data taken 
at LEP2 extends towards lower y cut values compared to the fit range for data taken at 
LEP1. This leads to an increased theoretical uncertainty for a s . Finally we combine all 



average y/s 
in GeV 


a s 


stat. 


exp. 


hadr. 


scale 


mass 


91 


0.1182 


0.0002 


0.0013 


0.0012 


0.0011 


0.0013 


133 


0.1067 


0.0042 


0.0045 


0.0010 


0.0011 


0.0010 


177 


0.1081 


0.0027 


0.0037 


0.0007 


0.0013 


0.0009 


197 


0.1070 


0.0013 


0.0039 


0.0006 


0.0016 


0.0009 



Table 2: The mean value of as for each energy interval, the statistical and experimental 
errors, and the errors due to hadronization and scale uncertainties and massless QCD. 



9 The uncertainty of the hadronization models and the detector correction using different Monte Carlo 
models is also less than 5%. 
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measurements to a single value of as(M z o), with the final result: 

a s (M z o) = 0.1182±0.0003(stat.)±0.0015(exp.)±0.0011(had.)±0.0012(scale)±0.0013(mass), 

(11) 

which is clearly dominated by the a s value obtained at 91 GeV. The weight as calculated 
in Eq. 10 for the data taken at 91 GeV is 84.3%. 

5.1.2 Scale Dependence of ct s 

For the fits in Section 5.1.1 the renormalization scale is set to the natural choice x^ = 
1. However, different schemes for the determination of the renormalization scale are 
proposed [32]. In this section we investigate two of these. 

In the optimized renormalization scheme, discussed in detail in [33], the minimization 
is performed with both as and treated as free parameters. However, at some low- 
statistics LEP2 energy points the fits did not converge. For this reason we applied a 
modified version of the optimized renormalization scheme. The optimized scale x° pt was 
determined with the high statistics sample at 91 GeV and then the same scale was used 
at higher energies. The result for rr° pt at 91 GeV is x° pt = 0.672 ± 0.033(stat.). The 
variation of x 2 /d.o.f. as a function of the scale x^ for 91 GeV data is shown in Figure 4. 
The combination of all energy points using the method described in Section 3.7 returns a 
value of 

a s (M z o) = 0.1187 ± 0.0005(stat.) ± 0.0018(exp.) ± 0.0011(had.) ± 0.0012(mass). (12) 

The weight as calculated in Eq. 10 for the data taken at 91 GeV is 73.1%, leading to an 
increased statistical uncertainty in the combined result. 

The second choice for the determination of the renormalization scale followed approx- 
imately the approach suggested by Stevenson [34]. The renormalization scale x™ m is 
specified by the point with a s having the least sensitivity to the renormalization scale 
x^. Using the 91 GeV data (Figure 4) x™ m is determined to be x° pt = 1.36. The fit is 
repeated at all other energy points with x^ set to x™" 1 . The combination of all energy 
points using the method described in Section 3.7 returns a value of 

a s (M z o) = 0.1182 ± 0.0003(stat.) ± 0.0014(exp.) ± 0.0012(had.) ± 0.0013(mass). (13) 

Both alternative choices for the renormalization scale return a value of as which is well 
within the variation of the systematic uncertainty due to missing higher order terms. 

5.2 Thrust minor and .D-parameter 

For the T m i n and D-parameter event-shape variables only the basic NLO calculations are 
available. For this reason and due to the problems discussed in the introduction we expect 
large systematic uncertainties due to missing higher order terms. The strong coupling a$ 
has not previously been determined using these four-jet event-shape variables. Here we 
perform a first measurement of as to quantify the various experimental, hadronization 
and theoretical uncertainties. 
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Studies with the basic NLO calculations were already performed with three-jet event- 
shape observables [33,35]. These investigations showed that fits with the renormalization 
scale fixed to x^ — 1 described the data worse than fits with both parameters as and 

free, which is usually interpreted to mean that the optimized scale-choice mimics the 
effect of higher order corrections. We use both approaches here. The D-parameter and 
T min event-shape distributions were published in [13] and are used here. For the QCD 
predictions we use the Debrecen 2.0 0(a$) calculations [14]. In the first approach the scale 
parameter x^, as discussed in Section 2.2, is set to one. The fit range is set from 0.010 
to 0.25 for the D-parameter and from 0.06 to 0.16 for T min . The fit range is determined 
by requiring that the hadronization and the detector correction be less than 50%. The 
results obtained at average centre-of-mass energies of 91 GeV and 197 GeV are shown in 
Figure 5. For the D-parameter the slope of the curve obtained from the theory prediction 
after the x 2 minimization does not give a good description of the data. Similar effects were 
already seen previously with three- and four-jet observables [8,10,33,35]. The results for 
the fit to the data between 91 and 209 GeV are combined to a single value for a s (M z o) 
using the method introduced in Section 3.7. The result for the D-parameter is 

a s (M z o) = 0.1047 ± 0.0014(stat.) ± 0.0037(exp.) ± 0.0061(had.) ± 0.0052(scale) (14) 

and for T min 

a s (M z o) = 0.1318 ± 0.0016(stat.) ± 0.0056(exp.) ± 0.0023(had.) ± O.Olll(scale), (15) 

both consistent with the result from the four-jet rate, albeit with large errors. All un- 
certainties are greater than those on the four-jet rate determination. In particular the 
systematic uncertainty due to higher order missing terms is significantly larger. 

The impact of the scale x^ on the result is studied with a revised fit using the mod- 
ified optimized renormalization scheme as in Section 5.1.2. The optimized scale x° pt is 
determined with 91 GeV data only and the fit repeated with ct s free and the scale set to 
x°p\ For the D-parameter the scale is set to x° pt = 0.0220 ± 0.0002(stat.) and for T min to 
rr° pt = 0.0089 ±0.0001(stat.) 10 . The description of the slope of the D-parameter improves 
significantly when the optimized renormalization scale parameter is used. The results are 
combined to a single value for a s (M z o) using the method introduced in Section 3.7. The 
result for the D-parameter is 

a s (M z o) = 0.0987 ± 0.0011(stat.) ± 0.0044(exp.) ± 0.0022(had.) (16) 

and for T min 

as (M z o) = 0.1100 ± 0.0011(stat.) ± 0.0045(exp.) ± 0.0032(had.). (17) 

Due to the problems discussed above we do not combine the as values based on the 
event-shape distributions with the one from the four-jet rate. The differences between 
the result evaluated at the natural scale and at the optimized scale suggests, as discussed 
in [7,8], that the uncalculated higher order terms are important. An improved theoret- 
ical understanding is necessary before as measurements using four-jet event-shapes are 
competitive. 

10 Other studies with observables using NLO predictions only with an optimized scale also result in 
values well below unity [10,33]. 
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6 Summary 



In this paper we present the measurements by the OPAL Collaboration of the strong 
coupling as from the four-jet rate at centre-of-mass energies between 91 and 209 GeV. 
The predictions of the PYTHIA, HERWIG and ARIADNE Monte Carlo models are found 
to be in general agreement with the measured distributions. 

From a fit of matched (9(o;g)+NLLA predictions to the four-jet rate, we have deter- 
mined the strong coupling ag(Mzo) = 0.1182 ± 0.0025(total error), which is consistent 
with the world average value of as(M z <>) = 0.1182 ± 0.0027 [36], based on a detailed 
evaluation of systematic errors. A fit to the D-parameter and the T min was performed 
using basic 0(a|) predictions. We find consistent results with significantly larger theo- 
retical uncertainties. Figure 6 summarizes the results obtained and compares them with 
results from previous OPAL publications using jet-rates [12], event shape distributions 
and moments [13]. The results are in good agreement with each other. Our result using 
the four-jet rate has the smallest uncertainty. This is mostly due to the comparatively 
small scale uncertainty, originating from the fact that the natural scale x M =l is close to 
the scale with the least sensitivity to as, as shown in Figure 4. 
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Figure 1: The four-jet rate distribution at hadron level as a function of the y cnt resolution 
parameter obtained with the Durham algorithm. The four-jet rates are shown for the 
data corrected to the hadron level at four average centre-of-mass energies between 91 and 
209 GeV together with predictions based on PYTHIA, HERWIG and ARIADNE Monte 
Carlo events generated at the averaged energy. The error bars show the statistical (inner 
error bars) and experimental uncertainties added in quadrature. Error bars not shown are 
smaller than the point size. The panel in each upper right corner shows the differences 
between data and Monte Carlo predictions, divided by the sum of the statistical and 
experimental error. For data points with no data events, the difference is set to zero. 
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Figure 2: The hadron-level four-jet rate distributions for energies between 91 GeV and 
209 GeV. The error bars show the statistical (inner error bars) and experimental uncer- 
tainties added in quadrature. When not shown, the errors are smaller than the point 
size. The curves show the theory prediction after \ 2 minimization within the fit range 
indicated. The fit range is determined as discussed in Section 5.1.1. The data points are 
strongly correlated and an enlarged fit range would lead to only a minor gain in statistical 
precision. 
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Figure 3: The values for as obtained by a fit to the four-jet rate with x M set to 1.0 in the 
four energy intervals. The error bars show the statistical (inner error bars) and the total 
error. The statistical error at 91 GeV is smaller than the point size. The lines indicate 
the current world average from [36] with the one standard deviation uncertainty. 
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Figure 4: The values of as and x 2 /d.o.f. from the fit to the four-jet rate at 91 GeV data 
as a function of the scale parameter x^. The arrows indicate the variation 0.5 < x^ < 2.0 
used to determine the theoretical systematic uncertainty. 
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Figure 5: The hadron-level /^-parameter and T m i n distributions for energies of 91 GeV 
and 197 GeV. The error bars show the statistical (inner error bars) and experimental 
uncertainties added in quadrature. When not shown, the errors are smaller than the 
point size. The curves indicate the theory prediction after ^-minimization with the 
renormalization scale set to 1.0. The fit range is determined as discussed in Section 5.2. 
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Figure 6: The results of the fit to the four-jet rate, the D-parameter and T m [ n . The 
value of the strong coupling as obtained from the four-jet rate represents our main result. 
The analyses of the D-parameter and T m i n return larger uncertainties as discussed in 
Section 5.2. Results published by the OPAL collaboration using different predictions and 
event shapes are shown as well. The value of as indicated by Jet Rates is obtained by 
[12]. The results of [13] are indicated as 'Event Shape Distributions' and 'Event Shape 
Moments'. The type of QCD prediction is indicated in the plot. The as world average is 
taken from [36]. 
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